Detailing the Regime Shifts in Maine Coastal Current Behavior
Published
August 21, 2025
Code
{library(sf) library(fvcom) library(tidyverse)library(gmRi)library(patchwork)library(rnaturalearth)library(showtext)library(ncdf4)# Cyclic color palettes in scico# From: https://www.fabiocrameri.ch/colourmaps/library(scico)library(legendry)library(ggh4x)library(ecodata)library(trend)}# namespace conflictsconflicted::conflict_prefer("select", "dplyr")conflicted::conflict_prefer("filter", "dplyr")# Set the themetheme_set(theme_gmri_simple() +theme(strip.text.y =element_text(angle =0),legend.position ="bottom", legend.title.position ="top", legend.title =element_text(hjust =0.5), strip.text =element_text(size =8),axis.text =element_text(size =8)))# Project pathslob_ecol_path <-cs_path("mills", "Projects/Lobster ECOL")fvcom_path <-cs_path("res", "FVCOM/Lobster-ECOL")poly_paths <-cs_path("mills", "Projects/Lobster ECOL/Spatial_Defs")# Shapefilesnew_england <-ne_states("united states of america", returnclass ="sf") %>%filter(postal %in%c("VT", "ME", "RI", "MA", "CT", "NH", "NY", "MD", "VA", "NJ", "DE", "NC", "PA", "WV"))canada <-ne_states("canada", returnclass ="sf")# # # Support functions for FVCOMsource(here::here("R/FVCOM_Support.R"))# New areas factor levelsareas_northsouth <-c("eastern maine", "central maine", "western maine", "eastern mass","southern mass", "rhode island shore","long island sound", "new jersey shore", "five fathom bank", "virginia shore","gom_gbk", "sne")
Code
# Use GMRI styleuse_gmri_style_rmd()
Code
# Stirnimann used these values in their paper:# l = 5, 10, 15, 17.5 years, with monthly data# Huber = 1# Subsampling = (l + 1) / 3# Load the function(s)source(here::here("rstars-master","rSTARS.R"))
Code
# Read Regions inproj_path <-cs_path("mills", "Projects/Lobster ECOL")# Load Shapefiles for inshore/offshorepoly_paths <-cs_path("mills", "Projects/Lobster ECOL/Spatial_Defs")# NEW Areas# clusters of statistical areas that align loosely with geography and management areasinshore_areas <-read_sf(str_c(poly_paths,"spatial_defs_2025/12nm_poly_statarea_merge.shp")) %>% janitor::clean_names() %>%mutate(area_type ="nearshore-coastal",area_id =tolower(short_name))# ecological production unitsoffshore_areas <-read_sf(str_c(poly_paths,"spatial_defs_2025/sne_gom_tocoast.shp")) %>% janitor::clean_names() %>%mutate(area_type ="offshore-regional",area_id =tolower(region))# Combine themstudy_regions <-bind_rows(st_transform(dplyr::select(inshore_areas, area_id, geometry), st_crs(offshore_areas)), dplyr::select(offshore_areas, area_id, geometry))
Reviewing STARS Regime Changes in Space/Time
This markdown reviews the various rstars regime shift results which were produced separately. I will begin at the largest geographic scales and work down to local timeseries:
Regime shifts for individual timeseries were tested using the STARS methodology. Any daily timeseries (temperature and salinity from ocean reanalysis models) were aggregated to a monthly temporal resolution, and any trends and seasonal cycles were removed.
The Marriott, Pope and Kendall (MPK) “pre-whitening” routine was used within the {rstars} algorithm to remove “red noise” (autoregressive processes, typically AR1) from the timeseries.
For more details on trend removal and pre-whitening methods see Rodionov 2006.
Shelf-Scale Ecodata Indices
There are 2-3 climate and oceanographic timeseries that operate at the broad regional scale of the Northeast shelf. These are:
The Gulf Stream Index (a metric indicating the North/South position of the Gulf Stream) based on SSH
The Northeast Channel Slopewater Proportions (the percentage of various water masses at the 150-200m depth entering GOM, using NERACOOS buoy data)
The North Atlantic Oscillation (atmospheric pressure differential between icelandic low and the Azores High)
These three indices affect conditions over large spatial scales and and are likely to either directly or indirectly impact more local scale changes.
These three metrics are available in the ecodata r package and can be pulled directly from the package.
Code
# GSI# Why is there more than one value per month?gsi <- ecodata::gsi %>%#glimpse()mutate(Time =as.Date(str_c(str_replace(Time, "[.]", "-"), "-01")))# use the old onegsi_old <- ecodata::gsi_old %>%#glimpse()mutate(Var ="gulf stream index old") %>%mutate(Time =as.Date(str_c(str_replace(Time, "[.]", "-"), "-01")))# NAOnao <- ecodata::nao%>%mutate(Time =as.Date(str_c(Time, "-01-01")))# Put them together to plotshelf_indices <-bind_rows(list(gsi, gsi_old, nao))# Plot the residuals from the trendggplot(shelf_indices, aes(Time, Value)) +geom_line(linewidth =0.6, alpha =0.8) +facet_grid(Var ~ ., scales ="free", labeller =label_wrap_gen(width=8)) +guides(color =guide_legend(nrow =2)) +labs(title ="Shelf Scale Ocean/Climate Metrics - Raw")
The Gulf Stream indices come as two monthly datasets, the other indices are annual.
Trends in Ecodata Indicators
A Mann-Kendall test can be used to determine whether there is any monotonic (increasing/decreasing) trends in a time series. I will be checking each timeseries along the way using this method to keep a tab on whether or not long-term trends existed, which may be relevant to ecosystem change regardless of stars regime results.
Code
# Perform test of monotonic trends and report as tableshelf_indices %>%split(.$Var) %>%map_dfr(function(x){# check pacing x_2000 <-filter(x, year(Time) ==2000) tempo <-if_else(nrow(x_2000) ==12, 12, 1)message(str_c(x$Var[[1]], " has ", tempo, " values")) trend <- trend::mk.test(x$Value)$p.value <0.05 rate <-ifelse( trend, coef(lm(Value ~ Time, data = x))[[2]] * tempo *10,NA)return(tibble("Trend?"= trend,"Decadal Rate"=round(rate, 3))) }, .id ="Var") %>% gt::gt() %>% gt::tab_header(title ="Ecodata Shelf-Scale Long Term Trends",subtitle ="Monotonic Trends Evaluated by Mann-Kendall Test") %>% gt::fmt_missing(columns =everything(), missing_text ="")
Ecodata Shelf-Scale Long Term Trends
Monotonic Trends Evaluated by Mann-Kendall Test
Var
Trend?
Decadal Rate
gulf stream index
TRUE
0.009
gulf stream index old
FALSE
north atlantic oscillation
FALSE
western gulf stream index
TRUE
0.008
The Gulf Stream and Western Gulf Stream indices show a long term increasing trend, with some evidence in the GSI that this quickly changed course around 2023.
Code
shelf_indices %>%filter(Var =="gulf stream index", year(Time) >2010) %>%ggplot(aes(Time, Value)) +geom_line() +geom_vline(xintercept =as.Date("2023-01-01"), linetype =3, color ="gray50") +facet_grid(Var~.) +labs(title ="Recent Course-Change in Gulf Stream Index",subtitle ="Rapid Gulf Stream Position Reset After 2023")
Because the presence of long-term trends can influence the results of tests for breakpoint/mean shifts, any long-term trends for each metric have been removed prior to regime shift tests on these metrics.
Code
# Run the shift test for summertime PCAshelf_indices_detrended <- shelf_indices %>%split(.$Var) %>%map_dfr(function(.x){# Detrend .x <- .x %>%arrange(Time) %>%mutate(time = Time,yr_num =year(time))# annual trend trend_mod <-lm(Value ~ Time, data = .x)# save the results .x <- broom::augment(x = trend_mod) %>%rename(trend_fit = .fitted,trend_resid = .resid) %>%full_join(.x, join_by(Time, Value)) %>%mutate(trend_resid =if_else(is.na(Value), NA, trend_resid))return(.x)})
The following plot shows what these timeseries look like after the removal of any monotonic trend:
Code
# Plot the residuals from the trendggplot(shelf_indices_detrended, aes(Time, trend_resid)) +geom_line(linewidth =0.6, alpha =0.8) +facet_grid(Var ~ ., scales ="free", labeller =label_wrap_gen(width=8)) +guides(color =guide_legend(nrow =2)) +labs(title ="Shelf Scale Ocean/Climate Metrics - Detrended")
Bacause the slopewater proportion contains NA values, we cannot evaluate it for breaks unless we impute missing values somehow or take a subset of time that is uninterrupted.
Code
# Run the regime shift testshelf_indices_rstars <- shelf_indices_detrended %>%split(.$Var) %>%map_dfr(function(.x){# Set the cutoff length, change it based on monthly/annual cutoff_length <-ifelse(str_detect(.x$Var[[1]], "index"),12*7,7)# This is only here because we have duplicate dates in the GSI .x <-distinct(.x, Time, .keep_all = T)# Get the results from that x_rstars <-rstars(data.timeseries =as.data.frame( .x[,c("Time", "trend_resid")]),l.cutoff = cutoff_length,pValue =0.05,Huber =1,Endfunction = T,preWhitening = T,OLS = F,MPK = T,IP4 = F,SubsampleSize = (cutoff_length +1)/3,returnResults = T) %>%mutate(EPU = .x$EPU[[1]],Value = .x$Value,shift_direction =case_when( RSI >0~"Shift Up", RSI <0~"Shift Down",TRUE~NA)) },.id ="Var" )
The results can be seen below:
Code
# Summarise the breakpoint locationsshelf_shift_points <- shelf_indices_rstars %>%filter(RSI !=0) %>% dplyr::select(Time, Var, EPU, shift_direction)# Plot the breaks over the monthly dataggplot() +# Add a vertical line using geom_segment with an arrowgeom_segment(data =filter( shelf_shift_points, str_detect(shift_direction, "Up")),linewidth =0.8,aes(x = Time, xend = Time, y =-Inf, , yend =Inf, color = shift_direction), arrow =arrow(length =unit(0.3, "cm"), ends ="last")) +geom_segment(data =filter( shelf_shift_points, !str_detect(shift_direction, "Up")),linewidth =0.8,aes(x = Time, xend = Time, y =-Inf, , yend =Inf, color = shift_direction), arrow =arrow(length =unit(0.3, "cm"), ends ="first")) +geom_line(data = shelf_indices_rstars,aes(Time, trend_resid),linewidth =0.4, alpha =0.5) +scale_color_gmri() +facet_grid(EPU * Var~., labeller =label_wrap_gen(width =8)) +scale_x_date(breaks =seq.Date(from =as.Date("1950-01-01"),to =as.Date("2020-01-01") ,by ="10 year"),limits =as.Date(c("1970-01-01", "2020-01-01")),labels = scales::label_date_short()) +theme(legend.position ="bottom",strip.text.y =element_text(angle =0)) +labs(color ="",y ="Measurement",x ="Date",title ="Shelf-Scale Ocean/Climate Metrics - STARS Changepoints",subtitle ="Performed on Detrended Timeseries")
Based on these results, there is evidence for breakpoints in the Gulf Stream Indices, and not in the NAO index.
For the EPU-scale metrics we have a number of indices available from the ecodata package:
The cold-pool index
The Northeast Channel Slopewater Proportion (from NERACOOS Buoy N)
Metrics of primary production and zooplankton community
Temperature and salinity timeseries specific to each area
Temperature and salinity is from either GLORYS or FVCOM, primary productivity is satellite derived (OC-CCI, SeaWiFS, MODIS-Aqua), and the zooplankton community indices are from the Gulf of Maine CPR transect.
Code
# # There are a ton here. We want primary productivity / chlor a, and maybe anomalies# chl_pp <- ecodata::chl_pp %>% # filter(str_detect(Var, "MONTHLY")) %>% # filter(str_detect(Var, "PPD|CHLOR_A")) %>% # separate(col = "Time", into = c("Period", "Time"), sep = "_") %>% # mutate(Time = as.Date(# str_c(# str_sub(Time, 1, 4),# str_sub(Time, 5, 6), # "01",# sep = "-")))# Annual will make life easierannual_chl_pp <- ecodata::annual_chl_pp %>%filter(str_detect(Var, "MEAN")) %>%separate(col ="Time", into =c("Period", "Time"), sep ="_") %>%mutate(Time =as.Date(str_c( Time,"01-01",sep ="-")))# Just take one cold pool index for nowcold_pool <- ecodata::cold_pool %>%#distinct(Source)filter(Var =="cold_pool_index") %>%mutate(Var =str_c(Source, Var, sep ="_"),Time =as.Date(str_c( Time,"01-01",sep ="-")))# Slopewaterslopewater <- ecodata::slopewater %>%mutate(Time =as.Date(str_c(Time, "-01-01"))) %>%filter(Time >as.Date("1990-01-01")) %>%drop_na()# Combine thoseepu_indices <-bind_rows(annual_chl_pp, slopewater, cold_pool) %>%group_by(Var, EPU) %>%arrange(Time) # Plot themggplot(epu_indices, aes(Time, Value)) +geom_line(linewidth =0.6, alpha =0.8) +facet_grid(Var * EPU ~ ., scales ="free", labeller =label_wrap_gen(width=8)) +guides(color =guide_legend(nrow =3)) +scale_x_date(breaks =seq.Date(from =as.Date("1950-01-01"),to =as.Date("2020-01-01") ,by ="10 year"),limits =as.Date(c("1970-01-01", "2020-01-01")),labels = scales::label_date_short()) +guides(color =guide_legend(nrow =2)) +labs(title ="EPU Scale Ocean/Ecological Metrics - Raw")
Long-Term EPU Scale Trends
Code
# Perform test of monotonic trends and report as tableepu_indices %>%mutate(var_epu =str_c(Var, EPU, sep ="X")) %>%split(.$var_epu) %>%map_dfr(function(x){# check pacing# Set the cutoff length, change it based on monthly/annual cutoff_length <-7# 7 years for annual trend <- trend::mk.test(x$Value)$p.value <0.05 rate <-ifelse( trend, coef(lm(Value ~ Time, data = x))[[2]] *12*10,NA)return(tibble("Trend?"= trend,"Decadal Rate"=round(rate, 3))) }, .id ="var_epu") %>%separate(var_epu, into =c("Var", "EPU"), sep ="X") %>% gt::gt() %>% gt::tab_header(title ="Ecodata Shelf-Scale Long Term Trends",subtitle ="Monotonic Trends Evaluated by Mann-Kendall Test") %>% gt::fmt_missing(columns =everything(), missing_text ="")
Ecodata Shelf-Scale Long Term Trends
Monotonic Trends Evaluated by Mann-Kendall Test
Var
EPU
Trend?
Decadal Rate
CHLOR_A_ANNUAL_MEAN
GB
TRUE
0.002
CHLOR_A_ANNUAL_MEAN
GOM
TRUE
0.004
CHLOR_A_ANNUAL_MEAN
MAB
FALSE
GLORYS_cold_pool_index
MAB
TRUE
0.016
LSLW proportion ne channel
GOM
FALSE
MOM6_cold_pool_index
MAB
FALSE
PPD_ANNUAL_MEAN
GB
TRUE
0.002
PPD_ANNUAL_MEAN
GOM
TRUE
0.002
PPD_ANNUAL_MEAN
MAB
TRUE
0.001
ROMS_cold_pool_index
MAB
TRUE
0.013
WSW proportion ne channel
GOM
FALSE
Detrending Ecodata Indicators
Most of these indicators have data at the monthly time-scale. To aid in regime change detection long-term year over year changes have been removed.
Once detrended, they can be checked for signs of regime changes.
Code
# Run rstars for those, temperature and salinity are done# Run the regime shift testepu_indices_rstars <- epu_indices_detrended %>%#filter(str_detect(Var, "proportion") == FALSE) %>% split(.$var_epu) %>%map_dfr(function(.x){# This is only here because we have duplicate dates in the GSI .x <-distinct(.x, Time, .keep_all = T)# cutoff length cutoff_length <-7# Get the results from that x_rstars <-rstars(data.timeseries =as.data.frame( .x[,c("Time", "trend_resid")]),l.cutoff = cutoff_length,pValue =0.05,Huber =1,Endfunction = T,preWhitening = T,OLS = F,MPK = T,IP4 = F,SubsampleSize = (cutoff_length +1)/3,returnResults = T) %>%mutate(Value = .x$Value,shift_direction =case_when( RSI >0~"Shift Up", RSI <0~"Shift Down",TRUE~NA)) },.id ="var_epu" ) %>%separate(var_epu, into =c("Var", "EPU"), sep ="X")
Primary Production and Cold-Pool Dynamics
The results from the STARS algorithm can be seen below:
Temperature and salinity timeseries from FVCOM were processed and had their STARS testing done separately in STARS_FVCOM.qmd. here are their results:
Code
# These are the raw monthly timeseries:# write_csv(surf_temp_monthly_shifts_raw, here::here("rstars_results/lobecol_stemp_monthly_shifts_raw.csv"))# write_csv(bot_temp_monthly_shifts_raw, here::here("rstars_results/lobecol_btemp_monthly_shifts_raw.csv"))# Load the data for the new regionsdaily_fvcom_temps <-read_csv(str_c(lob_ecol_path, "FVCOM_processed/area_timeseries/new_regions_fvcom_temperatures_daily.csv")) %>%mutate(yday = lubridate::yday(time),year = lubridate::year(time),month = lubridate::month(time),depth_type =if_else(area_id %in%c("gom_gbk", "sne"), "offshore", "nearshore"),area_id =factor(area_id, levels = areas_northsouth))# Run the monthly versionsmonthly_fvcom_temps <- daily_fvcom_temps %>%group_by(area_id, year, month, depth_type) %>%summarise(surface_t =mean(surface_t, na.rm = T),bottom_t =mean(bottom_t, na.rm = T),.groups ="drop") %>%mutate(yr_num =as.numeric(as.character(year)),month =factor(month),time =as.Date(str_c( year,str_pad(month, side ="left", pad ="0", width =2),"15", sep ="-"))) %>%rename(surface_temperature = surface_t,bottom_temperature = bottom_t) %>%pivot_longer(ends_with("temperature"), names_to ="var", values_to ="values")# Monthly Salinitymonthly_fvcom_sal <-read_csv(str_c(lob_ecol_path, "FVCOM_processed/area_timeseries/new_regions_fvcom_salinity_monthly_gom3.csv")) %>%mutate(yr_num =as.numeric(as.character(year(time))),month =factor(lubridate::month(time)),depth_type =if_else(area_id %in%c("gom_gbk", "sne"), "offshore", "nearshore"),area_id =factor(area_id, levels = areas_northsouth) ) %>%pivot_longer(ends_with("salinity"), names_to ="var", values_to ="values")# Combine thesefvcom_tempsal_raw <-bind_rows( monthly_fvcom_temps, monthly_fvcom_sal)
Temp/Sal Trends
Code
# Perform test of monotonic trends and report as tablefvcom_trends <- fvcom_tempsal_raw %>%mutate(var_area =str_c(var, area_id, sep ="X")) %>%split(.$var_area) %>%map_dfr(function(x){# check pacing# Set the cutoff length, change it based on monthly/annual cutoff_length <-7# 7 years for annual trend <- trend::mk.test(x$values)$p.value <0.05 rate <-ifelse( trend, coef(lm(values ~ time, data = x))[[2]] *12*10,NA)return(tibble("Trend?"= trend,"Decadal Rate"=round(rate, 3))) }, .id ="var_area") %>%separate(var_area, into =c("var", "area_id"), sep ="X") fvcom_trends %>%filter(#str_detect(var, "temperature"),`Trend?`, area_id %in%c("gom_gbk", "sne")) %>%group_by(area_id) %>% gt::gt() %>% gt::tab_header(title ="FVCOM Offshore Long Term Temperature Trends",subtitle ="Monotonic Trends Evaluated by Mann-Kendall Test") %>% gt::fmt_missing(columns =everything(), missing_text ="")
FVCOM Offshore Long Term Temperature Trends
Monotonic Trends Evaluated by Mann-Kendall Test
var
Trend?
Decadal Rate
sne
bottom_salinity
TRUE
0.002
bottom_temperature
TRUE
0.010
surface_temperature
TRUE
0.014
gom_gbk
bottom_temperature
TRUE
0.009
surface_salinity
TRUE
-0.001
surface_temperature
TRUE
0.013
Temp/Sal STARS Breaks
In STARS_FVCOM.qmd I used temperature timeseries to explore the impacts of performing regime shift testing on raw or detrended monthly data. These tests helped reinforce that the suggested preprocessing (detrending etc.) did help isolate step-changes in the timeseries which may be related to a regime change.
The following breaks were identified from these detrended series:
Code
# Load the regime shift results from STARS_FVCOM.qmd# These are from detrended datasurf_sal_monthly_shifts <-read_csv( here::here("rstars_results/lobecol_ssal_monthly_shifts_detrended.csv")) %>%mutate(Var ="Surface Salinity") %>%rename(detrended_vals = ssal_model_resid) bot_sal_monthly_shifts <-read_csv( here::here("rstars_results/lobecol_bsal_monthly_shifts_detrended.csv")) %>%mutate(Var ="Bottom Salinity") %>%rename(detrended_vals = bsal_model_resid)surf_temp_monthly_shifts <-read_csv( here::here("rstars_results/lobecol_stemp_monthly_shifts_detrended.csv")) %>%mutate(Var ="Surface Temperature") %>%rename(detrended_vals = stemp_model_resid)bot_temp_monthly_shifts <-read_csv( here::here("rstars_results/lobecol_btemp_monthly_shifts_detrended.csv")) %>%mutate(Var ="Bottom Temperature") %>%rename(detrended_vals = btemp_model_resid)# Put them togethertempsal <-bind_rows(list(surf_sal_monthly_shifts, bot_sal_monthly_shifts, surf_temp_monthly_shifts, bot_temp_monthly_shifts)) %>%mutate(shift_direction =case_when( RSI >0~"Shift Up", RSI <0~"Shift Down",TRUE~NA),area_id =factor(area_id, levels = areas_northsouth))offshore_tempsal <- tempsal %>%filter(area_id %in%tolower(c("GOM_GBK", "SNE")))inshore_tempsal <- tempsal %>%filter(area_id %in%tolower(c("GOM_GBK", "SNE")) ==FALSE)# Pull the shift pointstempsal_shifts <- tempsal %>%filter(RSI !=0)# Split surface and bottomoffshore_tempsal_shifts <- offshore_tempsal %>%filter(RSI !=0) %>% dplyr::select(time, Var, area_id, shift_direction)inshore_tempsal_shifts <- inshore_tempsal %>%filter(RSI !=0) %>% dplyr::select(time, Var, area_id, shift_direction)
Code
# Plot the offshore shifts# Plot the breaks over the monthly dataggplot() +# Add a vertical line using geom_segment with an arrowgeom_segment(data =filter( offshore_tempsal_shifts, str_detect(shift_direction, "Up")),linewidth =0.8,aes(x = time, xend = time, y =-Inf, yend =Inf, color = shift_direction), arrow =arrow(length =unit(0.3, "cm"), ends ="last")) +geom_segment(data =filter( offshore_tempsal_shifts, !str_detect(shift_direction, "Up")),linewidth =0.8,aes(x = time, xend = time, y =-Inf, yend =Inf, color = shift_direction), arrow =arrow(length =unit(0.3, "cm"), ends ="first")) +geom_line(data = offshore_tempsal,aes(time, detrended_vals),linewidth =0.4, alpha =0.5) +scale_color_gmri() +facet_nested(area_id * Var~., labeller =label_wrap_gen(width =8), scales ="free") +scale_x_date(breaks =seq.Date(from =as.Date("1950-01-01"),to =as.Date("2020-01-01") ,by ="10 year"),limits =as.Date(c("1978-01-01", "2020-01-01")),labels = scales::label_date_short()) +theme(legend.position ="bottom",strip.text.y =element_text(angle =0)) +labs(color ="",y ="Measurement",x ="Date",title ="FVCOM EPU Scale Temp/Sal Metrics",subtitle ="STARS Regime Shifts")
A change in SNE salinity appears to have occured around 1992.
Surface temperatures fell in SNE around 2002, but they rose again in 2011 along with GOM+GBK the same year.
Work by Andy Pershing helped develop an understanding that the Gulf of Mane’s zooplankton community in a given year is often one of two groups with different life history and size characteristics. There is a large copepod community, of which Calanus finmarchicus (a large bodied, lipid rich species) is prominent, and a second community which is composed of smaller-bodied and more opportunistic zooplankton species. These two communities compete for the same prey resources, and are typically out of phase with one-another. A principal component analysis using the continuous plankton recorded data has been used as a proxy for which community is dominant each year.
Code
# From Pershing & Kemberling# PC1 explains 53.62% of variance# PC2 explains 27.9%# PC1 is associated with centropages, oithona, para-pseudocalanus# PC2 is C. Fin, Metridia, & Euphausiacea# Load the CPR datacpr_community <-read_csv(here::here("local_data", "cpr_focal_pca_timeseries_period_1961-2017.csv")) %>%rename(PC1_small_zoo =`First Mode`,PC2_large_zoo =`Second Mode`) %>%select(-c(pca_period, taxa_used)) %>%pivot_longer(cols =starts_with("PC"), names_to ="Var", values_to ="value")
Taking PCA timeseries as proxies for those communities and evaluating them for breakpoints gives the following results.
Code
# Run breakpoints in CPR PCA# Run the regime shift testcpr_indices_rstars <- cpr_community %>%filter(year >1976) %>%mutate(EPU ="GOM",var_epu =str_c(Var, "X", EPU)) %>%split(.$var_epu) %>%map_dfr(function(.x){# detrend trend_mod <-lm(value ~ year, data = .x) .x$trend_resid <-resid(trend_mod)# cutoff length cutoff_length <-7# Get the results from that x_rstars <-rstars(data.timeseries =as.data.frame( .x[,c("year", "trend_resid")]),l.cutoff = cutoff_length,pValue =0.05,Huber =1,Endfunction = T,preWhitening = T,OLS = F,MPK = T,IP4 = F,SubsampleSize = (cutoff_length +1)/3,returnResults = T) %>%mutate(Value = .x$value,shift_direction =case_when( RSI >0~"Shift Up", RSI <0~"Shift Down",TRUE~NA)) },.id ="var_epu" ) %>%separate(var_epu, into =c("Var", "EPU"), sep ="X") %>%mutate(time =as.Date(str_c(year, "-01-01")))
# Abundance per 100m3 for different taxa# ecodata::zoo_regime %>% distinct(Var) %>% pull() %>% sort()ecomon_zoo <- ecodata::zoo_regimeecomon_zoo %>%filter(!str_detect(Var, "fish|clauso|gas")) %>%ggplot(aes(Time, Value)) +geom_line() +facet_grid(Var ~ EPU)
Code
# We might be able to just pull out the seven focal species and then repeat the PCA process...# Issues:# ecomon doesn't split the calanus into adult/juvenile# Para and Pseaudocalana are split# Euph also has a Euph1
Abundance anomalies are computed from the expected abundance on the day of sample collection. Abundance anomaly time series are constructed for Centropages typicus, Pseudocalanus spp., Calanus finmarchicus, and total zooplankton biovolume. The small-large copepod size index is computed by averaging the individual abundance anomalies of Pseudocalanus spp., Centropages hamatus, Centropages typicus, and Temora longicornis, and subtracting the abundance anomaly of Calanus finmarchicus. This index tracks the overall dominance of the small bodied copepods relative to the largest copepod in the Northeast U.S. region, Calanus finmarchicus.
Code
# This has "LgCopepods" & "SmCopepods", which could produce large/small index# ecodata::zoo_abundance_anom %>% distinct(Var) %>% pull() %>% sort()zoo_lg_small <- ecodata::zoo_abundance_anom %>%filter(Var %in%c("LgCopepods", "SmCopepods")) %>%mutate(Value =as.numeric(Value)) %>%pivot_wider(values_from ="Value", names_from ="Var") %>%mutate(small_large_index = SmCopepods - LgCopepods)ggplot(zoo_lg_small, aes(Time, small_large_index)) +geom_line() +geom_hline(yintercept =0) +facet_grid(EPU~., scales ="free") +labs(y ="Small-Large Copepod Index\n(More Large Copepods <-----> More Small Copepods)")
Code
# # This is too many vars# ecodata::zooplankton_index %>% distinct(Var) %>% pull()# # # Not informative mechanistically# ecodata::zoo_diversity# # # BOLD move on absolute abundances# ecodata::zoo_strat_abun
NEEDS: MCC & Lobster Predator Indices
There are two EPU-Scale indices that we need to develop. This is the MCC index, and a lobster predator abundance index.
The Gulf of Maine Coastal Current plays an important role in transporting lobster larva and their recruitment form year-to-year. The degree of “connected-ness” of the Western and Eastern portions of this current have been used in the past to inform expectations of lobster recruitment.
Local/Nearshore Shifts
Conditions closer to the coast show the following long-term trends:
Code
# Temperaturefvcom_trends %>%filter(str_detect(var, "temperature"),`Trend?`,!area_id %in%c("gom_gbk", "sne")) %>%group_by(area_id) %>% gt::gt() %>% gt::tab_header(title ="FVCOM Offshore Long Term Salinity Trends",subtitle ="Monotonic Trends Evaluated by Mann-Kendall Test") %>% gt::fmt_missing(columns =everything(), missing_text ="")
FVCOM Offshore Long Term Salinity Trends
Monotonic Trends Evaluated by Mann-Kendall Test
var
Trend?
Decadal Rate
eastern maine
bottom_temperature
TRUE
0.007
surface_temperature
TRUE
0.009
new jersey shore
bottom_temperature
TRUE
0.015
surface_temperature
TRUE
0.014
central maine
surface_temperature
TRUE
0.011
eastern mass
surface_temperature
TRUE
0.018
western maine
surface_temperature
TRUE
0.016
Code
# Salinityfvcom_trends %>%filter(!str_detect(var, "temperature"),`Trend?`,!area_id %in%c("gom_gbk", "sne")) %>%group_by(area_id) %>% gt::gt() %>% gt::tab_header(title ="FVCOM Offshore Long Term Salinity Trends",subtitle ="Monotonic Trends Evaluated by Mann-Kendall Test") %>% gt::fmt_missing(columns =everything(), missing_text ="")
FVCOM Offshore Long Term Salinity Trends
Monotonic Trends Evaluated by Mann-Kendall Test
var
Trend?
Decadal Rate
central maine
bottom_salinity
TRUE
-0.001
surface_salinity
TRUE
-0.003
eastern mass
bottom_salinity
TRUE
-0.004
surface_salinity
TRUE
-0.006
five fathom bank
bottom_salinity
TRUE
0.007
surface_salinity
TRUE
0.007
long island sound
bottom_salinity
TRUE
0.005
surface_salinity
TRUE
0.006
new jersey shore
bottom_salinity
TRUE
0.004
surface_salinity
TRUE
0.004
southern mass
bottom_salinity
TRUE
-0.003
surface_salinity
TRUE
-0.003
virginia shore
bottom_salinity
TRUE
0.003
surface_salinity
TRUE
0.003
western maine
bottom_salinity
TRUE
-0.002
surface_salinity
TRUE
-0.004
Temperature and Salinity Breaks
Code
# Plot the breaks over the monthly dataggplot() +# Add a vertical line using geom_segment with an arrowgeom_segment(data =filter( inshore_tempsal_shifts, str_detect(shift_direction, "Up")),linewidth =0.8,aes(x = time, xend = time, y =-Inf, yend =Inf, color = shift_direction), arrow =arrow(length =unit(0.3, "cm"), ends ="last")) +geom_segment(data =filter( inshore_tempsal_shifts, !str_detect(shift_direction, "Up")),linewidth =0.8,aes(x = time, xend = time, y =-Inf, yend =Inf, color = shift_direction), arrow =arrow(length =unit(0.3, "cm"), ends ="first")) +geom_line(data =filter(inshore_tempsal),aes(time, detrended_vals),linewidth =0.4, alpha =0.5) +scale_color_gmri() +facet_grid(area_id~Var, labeller =label_wrap_gen(width =8), scales ="free") +scale_x_date(breaks =seq.Date(from =as.Date("1950-01-01"),to =as.Date("2020-01-01") ,by ="10 year"),limits =as.Date(c("1978-01-01", "2020-01-01")),labels = scales::label_date_short()) +theme(legend.position ="bottom",strip.text.y =element_text(angle =0)) +labs(color ="",y ="Measurement",x ="Date",title ="FVCOM Inshore Scale Sal Metrics",subtitle ="STARS Regime Shifts")
Salinity
Code
# Plot the breaks over the monthly dataggplot() +# Add a vertical line using geom_segment with an arrowgeom_segment(data =filter( inshore_tempsal_shifts,str_detect(Var, "Salinity"),str_detect(shift_direction, "Up")),linewidth =0.8,aes(x = time, xend = time, y =-Inf, yend =Inf, color = shift_direction), arrow =arrow(length =unit(0.3, "cm"), ends ="last")) +geom_segment(data =filter( inshore_tempsal_shifts, str_detect(Var, "Salinity"),!str_detect(shift_direction, "Up")),linewidth =0.8,aes(x = time, xend = time, y =-Inf, yend =Inf, color = shift_direction), arrow =arrow(length =unit(0.3, "cm"), ends ="first")) +geom_line(data =filter(inshore_tempsal, str_detect(Var, "Salinity")),aes(time, detrended_vals),linewidth =0.4, alpha =0.5) +scale_color_gmri() +facet_grid(area_id~Var, labeller =label_wrap_gen(width =8), scales ="free") +scale_x_date(breaks =seq.Date(from =as.Date("1950-01-01"),to =as.Date("2020-01-01") ,by ="10 year"),limits =as.Date(c("1978-01-01", "2020-01-01")),labels = scales::label_date_short()) +theme(legend.position ="bottom",strip.text.y =element_text(angle =0)) +labs(color ="",y ="Salinity Anomaly",x ="Date",title ="FVCOM Hindcast Inshore Salinity",subtitle ="STARS Regime Shifts")
In addition to breaks in absolute temperatures, there is interest in the amount of time spent in favorable (12-18C) and unfavorable conditions (20C).
These use daily bottom temperatures:
Code
# Load the data for the new regionsdaily_fvcom_temps <-read_csv(str_c(lob_ecol_path, "FVCOM_processed/area_timeseries/new_regions_fvcom_temperatures_daily.csv")) %>%mutate(depth_type =if_else(area_id %in%c("gom_gbk", "SNE"), "offshore", "nearshore"),area_id =factor(area_id, levels = areas_northsouth) )# # Monthly Salinity# monthly_fvcom_sal <- read_csv(# str_c(lob_ecol_path, "FVCOM_processed/area_timeseries/new_regions_fvcom_salinity_monthly_gom3.csv")) %>% # mutate(# depth_type = if_else(area_id %in% c("gom_gbk", "SNE"), "offshore", "nearshore"),# area_id = factor(area_id, levels = areas_northsouth)# )# The degday functions were built with the ability to model daily cycles# the sine and triangle methods accomodate this# since we only have daily data the simple average is probably the way to do itthresh_low <-10thresh_up <-18# library(degday)# Get Monthly Totals in Rangesdd_monthly <- daily_fvcom_temps %>%mutate(year = lubridate::year(time),month = lubridate::month(time),opt_btemp =if_else(between(bottom_t, 10,18), 1, 0),stress_btemp =if_else(bottom_t >18, 1, NA),cold_btemp =if_else(bottom_t <10, 1, NA)) %>%group_by(area_id, year, month) %>%summarise(across(ends_with("temp"), ~sum(.x, na.rm = T)),.groups ="drop") %>%pivot_longer(cols =ends_with("temp"), names_to ="var", values_to ="totals") %>%mutate(time =as.Date(str_c(year,"-01-01")) +months(month-1),area_id =factor(area_id, levels = areas_northsouth),var =case_match( var,"opt_btemp"~"Preferred Bottom Temperatures 10-18C", "stress_btemp"~"Heat Stress Conditions >18C","cold_btemp"~"Below Preferred Conditions <10C") )
Code
# Do annual, Monthly values looked insanedd_annual <- dd_monthly %>%group_by(year, area_id, var) %>%summarise(across(totals, sum)) %>%mutate(var_area =str_c(var, area_id, sep ="X"),time =as.Date(str_c(year, "-01-01")))# Plot themdd_annual %>%ggplot() +geom_area(aes(year, y = totals, fill = var)) +scale_fill_manual(values =c("lightblue", "#ea4f12", "#057872")) +facet_grid(area_id~var) +scale_x_continuous(expand =expansion(add =c(0,0))) +theme(strip.text.y =element_text(angle =0),legend.position ="bottom") +guides(fill =guide_legend(nrow =2,title.position ="top",title.hjust =0.5))+labs(y ="Days in Range",fill ="Daily Temperature Conditions", color ="",title ="FVCOM Bottom Temperature Degree-Days")
Summary table should have the region, the variable, whether there is a trend or not, what the rate is, then whether there have been breakpoints, and when they were (mm/yyyy).
Source Code
---title: "Regime Shift Review"description: | Detailing the Regime Shifts in Maine Coastal Current Behaviordate: "Updated on: `r Sys.Date()`"format: html: code-fold: true code-tools: true df-print: kable self-contained: trueexecute: echo: true warning: false message: false fig.align: center comment: ""---```{r}{library(sf) library(fvcom) library(tidyverse)library(gmRi)library(patchwork)library(rnaturalearth)library(showtext)library(ncdf4)# Cyclic color palettes in scico# From: https://www.fabiocrameri.ch/colourmaps/library(scico)library(legendry)library(ggh4x)library(ecodata)library(trend)}# namespace conflictsconflicted::conflict_prefer("select", "dplyr")conflicted::conflict_prefer("filter", "dplyr")# Set the themetheme_set(theme_gmri_simple() +theme(strip.text.y =element_text(angle =0),legend.position ="bottom", legend.title.position ="top", legend.title =element_text(hjust =0.5), strip.text =element_text(size =8),axis.text =element_text(size =8)))# Project pathslob_ecol_path <-cs_path("mills", "Projects/Lobster ECOL")fvcom_path <-cs_path("res", "FVCOM/Lobster-ECOL")poly_paths <-cs_path("mills", "Projects/Lobster ECOL/Spatial_Defs")# Shapefilesnew_england <-ne_states("united states of america", returnclass ="sf") %>%filter(postal %in%c("VT", "ME", "RI", "MA", "CT", "NH", "NY", "MD", "VA", "NJ", "DE", "NC", "PA", "WV"))canada <-ne_states("canada", returnclass ="sf")# # # Support functions for FVCOMsource(here::here("R/FVCOM_Support.R"))# New areas factor levelsareas_northsouth <-c("eastern maine", "central maine", "western maine", "eastern mass","southern mass", "rhode island shore","long island sound", "new jersey shore", "five fathom bank", "virginia shore","gom_gbk", "sne")``````{r}#| label: style-sheet#| results: asis# Use GMRI styleuse_gmri_style_rmd()``````{r}#| label: load rstars functions# Stirnimann used these values in their paper:# l = 5, 10, 15, 17.5 years, with monthly data# Huber = 1# Subsampling = (l + 1) / 3# Load the function(s)source(here::here("rstars-master","rSTARS.R"))``````{r}#| label: load shapes# Read Regions inproj_path <-cs_path("mills", "Projects/Lobster ECOL")# Load Shapefiles for inshore/offshorepoly_paths <-cs_path("mills", "Projects/Lobster ECOL/Spatial_Defs")# NEW Areas# clusters of statistical areas that align loosely with geography and management areasinshore_areas <-read_sf(str_c(poly_paths,"spatial_defs_2025/12nm_poly_statarea_merge.shp")) %>% janitor::clean_names() %>%mutate(area_type ="nearshore-coastal",area_id =tolower(short_name))# ecological production unitsoffshore_areas <-read_sf(str_c(poly_paths,"spatial_defs_2025/sne_gom_tocoast.shp")) %>% janitor::clean_names() %>%mutate(area_type ="offshore-regional",area_id =tolower(region))# Combine themstudy_regions <-bind_rows(st_transform(dplyr::select(inshore_areas, area_id, geometry), st_crs(offshore_areas)), dplyr::select(offshore_areas, area_id, geometry))``````{r}#| label: fonts-config#| echo: false# Path to the directory containing the font file (replace with your actual path)font_dir <-paste0(system.file("stylesheets", package ="gmRi"), "/GMRI_fonts/Avenir/")# Register the fontfont_add(family ="Avenir",file.path(font_dir, "LTe50342.ttf"),bold =file.path(font_dir, "LTe50340.ttf"),italic =file.path(font_dir, "LTe50343.ttf"),bolditalic =file.path(font_dir, "LTe50347.ttf"))# Load the fontshowtext::showtext_auto()```# Reviewing STARS Regime Changes in Space/TimeThis markdown reviews the various rstars regime shift results which were produced separately. I will begin at the largest geographic scales and work down to local timeseries:Regime shifts for individual timeseries were tested using the STARS methodology. Any daily timeseries (temperature and salinity from ocean reanalysis models) were aggregated to a monthly temporal resolution, and any trends and seasonal cycles were removed.The Marriott, Pope and Kendall (MPK) "pre-whitening" routine was used within the {rstars} algorithm to remove "red noise" (autoregressive processes, typically AR1) from the timeseries.For more details on trend removal and pre-whitening methods see Rodionov 2006.# Shelf-Scale Ecodata IndicesThere are 2-3 climate and oceanographic timeseries that operate at the broad regional scale of the Northeast shelf. These are:1. The Gulf Stream Index (a metric indicating the North/South position of the Gulf Stream) based on SSH2. The Northeast Channel Slopewater Proportions (the percentage of various water masses at the 150-200m depth entering GOM, using NERACOOS buoy data)3. The North Atlantic Oscillation (atmospheric pressure differential between icelandic low and the Azores High)These three indices affect conditions over large spatial scales and and are likely to either directly or indirectly impact more local scale changes.These three metrics are available in the `ecodata` r package and can be pulled directly from the package.```{r}#| label: load shef-scale indices# GSI# Why is there more than one value per month?gsi <- ecodata::gsi %>%#glimpse()mutate(Time =as.Date(str_c(str_replace(Time, "[.]", "-"), "-01")))# use the old onegsi_old <- ecodata::gsi_old %>%#glimpse()mutate(Var ="gulf stream index old") %>%mutate(Time =as.Date(str_c(str_replace(Time, "[.]", "-"), "-01")))# NAOnao <- ecodata::nao%>%mutate(Time =as.Date(str_c(Time, "-01-01")))# Put them together to plotshelf_indices <-bind_rows(list(gsi, gsi_old, nao))# Plot the residuals from the trendggplot(shelf_indices, aes(Time, Value)) +geom_line(linewidth =0.6, alpha =0.8) +facet_grid(Var ~ ., scales ="free", labeller =label_wrap_gen(width=8)) +guides(color =guide_legend(nrow =2)) +labs(title ="Shelf Scale Ocean/Climate Metrics - Raw")```The Gulf Stream indices come as two monthly datasets, the other indices are annual. ### Trends in Ecodata IndicatorsA Mann-Kendall test can be used to determine whether there is any monotonic (increasing/decreasing) trends in a time series. I will be checking each timeseries along the way using this method to keep a tab on whether or not long-term trends existed, which may be relevant to ecosystem change regardless of stars regime results.```{r}# Perform test of monotonic trends and report as tableshelf_indices %>%split(.$Var) %>%map_dfr(function(x){# check pacing x_2000 <-filter(x, year(Time) ==2000) tempo <-if_else(nrow(x_2000) ==12, 12, 1)message(str_c(x$Var[[1]], " has ", tempo, " values")) trend <- trend::mk.test(x$Value)$p.value <0.05 rate <-ifelse( trend, coef(lm(Value ~ Time, data = x))[[2]] * tempo *10,NA)return(tibble("Trend?"= trend,"Decadal Rate"=round(rate, 3))) }, .id ="Var") %>% gt::gt() %>% gt::tab_header(title ="Ecodata Shelf-Scale Long Term Trends",subtitle ="Monotonic Trends Evaluated by Mann-Kendall Test") %>% gt::fmt_missing(columns =everything(), missing_text ="")```The Gulf Stream and Western Gulf Stream indices show a long term increasing trend, with some evidence in the GSI that this quickly changed course around 2023.```{r}shelf_indices %>%filter(Var =="gulf stream index", year(Time) >2010) %>%ggplot(aes(Time, Value)) +geom_line() +geom_vline(xintercept =as.Date("2023-01-01"), linetype =3, color ="gray50") +facet_grid(Var~.) +labs(title ="Recent Course-Change in Gulf Stream Index",subtitle ="Rapid Gulf Stream Position Reset After 2023")```Because the presence of long-term trends can influence the results of tests for breakpoint/mean shifts, any long-term trends for each metric have been removed prior to regime shift tests on these metrics.```{r}#| fig-height: 5# Run the shift test for summertime PCAshelf_indices_detrended <- shelf_indices %>%split(.$Var) %>%map_dfr(function(.x){# Detrend .x <- .x %>%arrange(Time) %>%mutate(time = Time,yr_num =year(time))# annual trend trend_mod <-lm(Value ~ Time, data = .x)# save the results .x <- broom::augment(x = trend_mod) %>%rename(trend_fit = .fitted,trend_resid = .resid) %>%full_join(.x, join_by(Time, Value)) %>%mutate(trend_resid =if_else(is.na(Value), NA, trend_resid))return(.x)})```The following plot shows what these timeseries look like after the removal of any monotonic trend:```{r}# Plot the residuals from the trendggplot(shelf_indices_detrended, aes(Time, trend_resid)) +geom_line(linewidth =0.6, alpha =0.8) +facet_grid(Var ~ ., scales ="free", labeller =label_wrap_gen(width=8)) +guides(color =guide_legend(nrow =2)) +labs(title ="Shelf Scale Ocean/Climate Metrics - Detrended")```Bacause the slopewater proportion contains NA values, we cannot evaluate it for breaks unless we impute missing values somehow or take a subset of time that is uninterrupted.```{r}# Run the regime shift testshelf_indices_rstars <- shelf_indices_detrended %>%split(.$Var) %>%map_dfr(function(.x){# Set the cutoff length, change it based on monthly/annual cutoff_length <-ifelse(str_detect(.x$Var[[1]], "index"),12*7,7)# This is only here because we have duplicate dates in the GSI .x <-distinct(.x, Time, .keep_all = T)# Get the results from that x_rstars <-rstars(data.timeseries =as.data.frame( .x[,c("Time", "trend_resid")]),l.cutoff = cutoff_length,pValue =0.05,Huber =1,Endfunction = T,preWhitening = T,OLS = F,MPK = T,IP4 = F,SubsampleSize = (cutoff_length +1)/3,returnResults = T) %>%mutate(EPU = .x$EPU[[1]],Value = .x$Value,shift_direction =case_when( RSI >0~"Shift Up", RSI <0~"Shift Down",TRUE~NA)) },.id ="Var" )```The results can be seen below:```{r}# Summarise the breakpoint locationsshelf_shift_points <- shelf_indices_rstars %>%filter(RSI !=0) %>% dplyr::select(Time, Var, EPU, shift_direction)# Plot the breaks over the monthly dataggplot() +# Add a vertical line using geom_segment with an arrowgeom_segment(data =filter( shelf_shift_points, str_detect(shift_direction, "Up")),linewidth =0.8,aes(x = Time, xend = Time, y =-Inf, , yend =Inf, color = shift_direction), arrow =arrow(length =unit(0.3, "cm"), ends ="last")) +geom_segment(data =filter( shelf_shift_points, !str_detect(shift_direction, "Up")),linewidth =0.8,aes(x = Time, xend = Time, y =-Inf, , yend =Inf, color = shift_direction), arrow =arrow(length =unit(0.3, "cm"), ends ="first")) +geom_line(data = shelf_indices_rstars,aes(Time, trend_resid),linewidth =0.4, alpha =0.5) +scale_color_gmri() +facet_grid(EPU * Var~., labeller =label_wrap_gen(width =8)) +scale_x_date(breaks =seq.Date(from =as.Date("1950-01-01"),to =as.Date("2020-01-01") ,by ="10 year"),limits =as.Date(c("1970-01-01", "2020-01-01")),labels = scales::label_date_short()) +theme(legend.position ="bottom",strip.text.y =element_text(angle =0)) +labs(color ="",y ="Measurement",x ="Date",title ="Shelf-Scale Ocean/Climate Metrics - STARS Changepoints",subtitle ="Performed on Detrended Timeseries")```Based on these results, there is evidence for breakpoints in the Gulf Stream Indices, and not in the NAO index.```{r}shelf_shift_points %>%group_by(Var) %>%arrange(Time, EPU) %>% gt::gt() %>% gt::tab_header(title ="Shelf-Scale Breaks")```# EPU-Scale Ecodata IndicesFor the EPU-scale metrics we have a number of indices available from the `ecodata` package:1. The cold-pool index2. The Northeast Channel Slopewater Proportion (from NERACOOS Buoy N)3. Metrics of primary production and zooplankton community4. Temperature and salinity timeseries specific to each areaTemperature and salinity is from either GLORYS or FVCOM, primary productivity is satellite derived (OC-CCI, SeaWiFS, MODIS-Aqua), and the zooplankton community indices are from the Gulf of Maine CPR transect.```{r}# # There are a ton here. We want primary productivity / chlor a, and maybe anomalies# chl_pp <- ecodata::chl_pp %>% # filter(str_detect(Var, "MONTHLY")) %>% # filter(str_detect(Var, "PPD|CHLOR_A")) %>% # separate(col = "Time", into = c("Period", "Time"), sep = "_") %>% # mutate(Time = as.Date(# str_c(# str_sub(Time, 1, 4),# str_sub(Time, 5, 6), # "01",# sep = "-")))# Annual will make life easierannual_chl_pp <- ecodata::annual_chl_pp %>%filter(str_detect(Var, "MEAN")) %>%separate(col ="Time", into =c("Period", "Time"), sep ="_") %>%mutate(Time =as.Date(str_c( Time,"01-01",sep ="-")))# Just take one cold pool index for nowcold_pool <- ecodata::cold_pool %>%#distinct(Source)filter(Var =="cold_pool_index") %>%mutate(Var =str_c(Source, Var, sep ="_"),Time =as.Date(str_c( Time,"01-01",sep ="-")))# Slopewaterslopewater <- ecodata::slopewater %>%mutate(Time =as.Date(str_c(Time, "-01-01"))) %>%filter(Time >as.Date("1990-01-01")) %>%drop_na()# Combine thoseepu_indices <-bind_rows(annual_chl_pp, slopewater, cold_pool) %>%group_by(Var, EPU) %>%arrange(Time) # Plot themggplot(epu_indices, aes(Time, Value)) +geom_line(linewidth =0.6, alpha =0.8) +facet_grid(Var * EPU ~ ., scales ="free", labeller =label_wrap_gen(width=8)) +guides(color =guide_legend(nrow =3)) +scale_x_date(breaks =seq.Date(from =as.Date("1950-01-01"),to =as.Date("2020-01-01") ,by ="10 year"),limits =as.Date(c("1970-01-01", "2020-01-01")),labels = scales::label_date_short()) +guides(color =guide_legend(nrow =2)) +labs(title ="EPU Scale Ocean/Ecological Metrics - Raw")```### Long-Term EPU Scale Trends```{r}# Perform test of monotonic trends and report as tableepu_indices %>%mutate(var_epu =str_c(Var, EPU, sep ="X")) %>%split(.$var_epu) %>%map_dfr(function(x){# check pacing# Set the cutoff length, change it based on monthly/annual cutoff_length <-7# 7 years for annual trend <- trend::mk.test(x$Value)$p.value <0.05 rate <-ifelse( trend, coef(lm(Value ~ Time, data = x))[[2]] *12*10,NA)return(tibble("Trend?"= trend,"Decadal Rate"=round(rate, 3))) }, .id ="var_epu") %>%separate(var_epu, into =c("Var", "EPU"), sep ="X") %>% gt::gt() %>% gt::tab_header(title ="Ecodata Shelf-Scale Long Term Trends",subtitle ="Monotonic Trends Evaluated by Mann-Kendall Test") %>% gt::fmt_missing(columns =everything(), missing_text ="")```### Detrending Ecodata IndicatorsMost of these indicators have data at the monthly time-scale. To aid in regime change detection long-term year over year changes have been removed.```{r}#| fig.height: 6# Detrend the epu stuffepu_indices_detrended <- epu_indices %>%mutate(var_epu =str_c(Var, EPU, sep ="X")) %>%split(.$var_epu) %>%map_dfr(function(.x){# Detrend .x <- .x %>%arrange(Time) %>%mutate(time = Time,yr_num =year(time))# annual trend trend_mod <-lm(Value ~ Time, data = .x)# save the results .x <- broom::augment(x = trend_mod) %>%rename(trend_fit = .fitted,trend_resid = .resid) %>%full_join(.x, join_by(Time, Value)) %>%mutate(trend_resid =if_else(is.na(Value), NA, trend_resid))return(.x)})# Plot detrendedggplot(epu_indices_detrended, aes(Time, trend_resid)) +geom_line(linewidth =0.6, alpha =0.8) +facet_grid(EPU * Var ~ ., scales ="free", labeller =label_wrap_gen(width=8)) +guides(color =guide_legend(nrow =3)) +scale_x_date(breaks =seq.Date(from =as.Date("1950-01-01"),to =as.Date("2020-01-01") ,by ="10 year"),limits =as.Date(c("1970-01-01", "2020-01-01")),labels = scales::label_date_short()) +labs(title ="EPU Scale Ocean/Climate Metrics - Detrended",y ="Metric")```Once detrended, they can be checked for signs of regime changes.```{r}# Run rstars for those, temperature and salinity are done# Run the regime shift testepu_indices_rstars <- epu_indices_detrended %>%#filter(str_detect(Var, "proportion") == FALSE) %>% split(.$var_epu) %>%map_dfr(function(.x){# This is only here because we have duplicate dates in the GSI .x <-distinct(.x, Time, .keep_all = T)# cutoff length cutoff_length <-7# Get the results from that x_rstars <-rstars(data.timeseries =as.data.frame( .x[,c("Time", "trend_resid")]),l.cutoff = cutoff_length,pValue =0.05,Huber =1,Endfunction = T,preWhitening = T,OLS = F,MPK = T,IP4 = F,SubsampleSize = (cutoff_length +1)/3,returnResults = T) %>%mutate(Value = .x$Value,shift_direction =case_when( RSI >0~"Shift Up", RSI <0~"Shift Down",TRUE~NA)) },.id ="var_epu" ) %>%separate(var_epu, into =c("Var", "EPU"), sep ="X")```### Primary Production and Cold-Pool DynamicsThe results from the STARS algorithm can be seen below:```{r}# Summarise the breakpoint locationsepu_shift_points <- epu_indices_rstars %>%filter(RSI !=0) %>% dplyr::select(Time, Var, EPU, shift_direction)# Plot the breaks over the monthly dataggplot() +geom_vline(data = epu_shift_points,aes(xintercept = Time,color = shift_direction),linewidth =1.5) +geom_line(data = epu_indices_rstars,aes(Time, Value),linewidth =0.4, alpha =0.5) +scale_color_gmri() +facet_grid(EPU * Var~., labeller =label_wrap_gen(width =8), scales ="free") +scale_x_date(breaks =seq.Date(from =as.Date("1950-01-01"),to =as.Date("2020-01-01") ,by ="10 year"),limits =as.Date(c("1970-01-01", "2020-01-01")),labels = scales::label_date_short()) +theme(legend.position ="bottom",strip.text.y =element_text(angle =0)) +labs(color ="",y ="Measurement",x ="Date",title ="EPU Scale Ocean/Climate Metrics - Detrended",subtitle ="STARS Regime Shifts")```Based on these results, we see no breakpoints in EPU-Scale measures of primary production or the cold pool dynamics.```{r}epu_shift_points %>%group_by(Var) %>%arrange(Time, EPU) %>% gt::gt() %>% gt::tab_header(title ="Ecodata EPU Scale Breakpoints")```### EPU-Scale Temperature and SalinityTemperature and salinity timeseries from FVCOM were processed and had their STARS testing done separately in `STARS_FVCOM.qmd`. here are their results:```{r}#| label: Load FVCOM Temperature and Salinity# These are the raw monthly timeseries:# write_csv(surf_temp_monthly_shifts_raw, here::here("rstars_results/lobecol_stemp_monthly_shifts_raw.csv"))# write_csv(bot_temp_monthly_shifts_raw, here::here("rstars_results/lobecol_btemp_monthly_shifts_raw.csv"))# Load the data for the new regionsdaily_fvcom_temps <-read_csv(str_c(lob_ecol_path, "FVCOM_processed/area_timeseries/new_regions_fvcom_temperatures_daily.csv")) %>%mutate(yday = lubridate::yday(time),year = lubridate::year(time),month = lubridate::month(time),depth_type =if_else(area_id %in%c("gom_gbk", "sne"), "offshore", "nearshore"),area_id =factor(area_id, levels = areas_northsouth))# Run the monthly versionsmonthly_fvcom_temps <- daily_fvcom_temps %>%group_by(area_id, year, month, depth_type) %>%summarise(surface_t =mean(surface_t, na.rm = T),bottom_t =mean(bottom_t, na.rm = T),.groups ="drop") %>%mutate(yr_num =as.numeric(as.character(year)),month =factor(month),time =as.Date(str_c( year,str_pad(month, side ="left", pad ="0", width =2),"15", sep ="-"))) %>%rename(surface_temperature = surface_t,bottom_temperature = bottom_t) %>%pivot_longer(ends_with("temperature"), names_to ="var", values_to ="values")# Monthly Salinitymonthly_fvcom_sal <-read_csv(str_c(lob_ecol_path, "FVCOM_processed/area_timeseries/new_regions_fvcom_salinity_monthly_gom3.csv")) %>%mutate(yr_num =as.numeric(as.character(year(time))),month =factor(lubridate::month(time)),depth_type =if_else(area_id %in%c("gom_gbk", "sne"), "offshore", "nearshore"),area_id =factor(area_id, levels = areas_northsouth) ) %>%pivot_longer(ends_with("salinity"), names_to ="var", values_to ="values")# Combine thesefvcom_tempsal_raw <-bind_rows( monthly_fvcom_temps, monthly_fvcom_sal)```#### Temp/Sal Trends```{r}# Perform test of monotonic trends and report as tablefvcom_trends <- fvcom_tempsal_raw %>%mutate(var_area =str_c(var, area_id, sep ="X")) %>%split(.$var_area) %>%map_dfr(function(x){# check pacing# Set the cutoff length, change it based on monthly/annual cutoff_length <-7# 7 years for annual trend <- trend::mk.test(x$values)$p.value <0.05 rate <-ifelse( trend, coef(lm(values ~ time, data = x))[[2]] *12*10,NA)return(tibble("Trend?"= trend,"Decadal Rate"=round(rate, 3))) }, .id ="var_area") %>%separate(var_area, into =c("var", "area_id"), sep ="X") fvcom_trends %>%filter(#str_detect(var, "temperature"),`Trend?`, area_id %in%c("gom_gbk", "sne")) %>%group_by(area_id) %>% gt::gt() %>% gt::tab_header(title ="FVCOM Offshore Long Term Temperature Trends",subtitle ="Monotonic Trends Evaluated by Mann-Kendall Test") %>% gt::fmt_missing(columns =everything(), missing_text ="")```#### Temp/Sal STARS BreaksIn `STARS_FVCOM.qmd` I used temperature timeseries to explore the impacts of performing regime shift testing on raw or detrended monthly data. These tests helped reinforce that the suggested preprocessing (detrending etc.) did help isolate step-changes in the timeseries which may be related to a regime change.The following breaks were identified from these detrended series:```{r}#| label: load detrended rsi for sal/temp# Load the regime shift results from STARS_FVCOM.qmd# These are from detrended datasurf_sal_monthly_shifts <-read_csv( here::here("rstars_results/lobecol_ssal_monthly_shifts_detrended.csv")) %>%mutate(Var ="Surface Salinity") %>%rename(detrended_vals = ssal_model_resid) bot_sal_monthly_shifts <-read_csv( here::here("rstars_results/lobecol_bsal_monthly_shifts_detrended.csv")) %>%mutate(Var ="Bottom Salinity") %>%rename(detrended_vals = bsal_model_resid)surf_temp_monthly_shifts <-read_csv( here::here("rstars_results/lobecol_stemp_monthly_shifts_detrended.csv")) %>%mutate(Var ="Surface Temperature") %>%rename(detrended_vals = stemp_model_resid)bot_temp_monthly_shifts <-read_csv( here::here("rstars_results/lobecol_btemp_monthly_shifts_detrended.csv")) %>%mutate(Var ="Bottom Temperature") %>%rename(detrended_vals = btemp_model_resid)# Put them togethertempsal <-bind_rows(list(surf_sal_monthly_shifts, bot_sal_monthly_shifts, surf_temp_monthly_shifts, bot_temp_monthly_shifts)) %>%mutate(shift_direction =case_when( RSI >0~"Shift Up", RSI <0~"Shift Down",TRUE~NA),area_id =factor(area_id, levels = areas_northsouth))offshore_tempsal <- tempsal %>%filter(area_id %in%tolower(c("GOM_GBK", "SNE")))inshore_tempsal <- tempsal %>%filter(area_id %in%tolower(c("GOM_GBK", "SNE")) ==FALSE)# Pull the shift pointstempsal_shifts <- tempsal %>%filter(RSI !=0)# Split surface and bottomoffshore_tempsal_shifts <- offshore_tempsal %>%filter(RSI !=0) %>% dplyr::select(time, Var, area_id, shift_direction)inshore_tempsal_shifts <- inshore_tempsal %>%filter(RSI !=0) %>% dplyr::select(time, Var, area_id, shift_direction)``````{r}#| fig.height: 5# Plot the offshore shifts# Plot the breaks over the monthly dataggplot() +# Add a vertical line using geom_segment with an arrowgeom_segment(data =filter( offshore_tempsal_shifts, str_detect(shift_direction, "Up")),linewidth =0.8,aes(x = time, xend = time, y =-Inf, yend =Inf, color = shift_direction), arrow =arrow(length =unit(0.3, "cm"), ends ="last")) +geom_segment(data =filter( offshore_tempsal_shifts, !str_detect(shift_direction, "Up")),linewidth =0.8,aes(x = time, xend = time, y =-Inf, yend =Inf, color = shift_direction), arrow =arrow(length =unit(0.3, "cm"), ends ="first")) +geom_line(data = offshore_tempsal,aes(time, detrended_vals),linewidth =0.4, alpha =0.5) +scale_color_gmri() +facet_nested(area_id * Var~., labeller =label_wrap_gen(width =8), scales ="free") +scale_x_date(breaks =seq.Date(from =as.Date("1950-01-01"),to =as.Date("2020-01-01") ,by ="10 year"),limits =as.Date(c("1978-01-01", "2020-01-01")),labels = scales::label_date_short()) +theme(legend.position ="bottom",strip.text.y =element_text(angle =0)) +labs(color ="",y ="Measurement",x ="Date",title ="FVCOM EPU Scale Temp/Sal Metrics",subtitle ="STARS Regime Shifts")```A change in SNE salinity appears to have occured around 1992.Surface temperatures fell in SNE around 2002, but they rose again in 2011 along with GOM+GBK the same year.```{r}offshore_tempsal_shifts %>%group_by(Var) %>%arrange(time, area_id) %>% gt::gt() %>% gt::tab_header(title ="EPU Scale Temp+Sal Breaks")```### CPR Community PCA IndexWork by Andy Pershing helped develop an understanding that the Gulf of Mane's zooplankton community in a given year is often one of two groups with different life history and size characteristics. There is a large copepod community, of which Calanus finmarchicus (a large bodied, lipid rich species) is prominent, and a second community which is composed of smaller-bodied and more opportunistic zooplankton species. These two communities compete for the same prey resources, and are typically out of phase with one-another. A principal component analysis using the continuous plankton recorded data has been used as a proxy for which community is dominant each year.```{r}# From Pershing & Kemberling# PC1 explains 53.62% of variance# PC2 explains 27.9%# PC1 is associated with centropages, oithona, para-pseudocalanus# PC2 is C. Fin, Metridia, & Euphausiacea# Load the CPR datacpr_community <-read_csv(here::here("local_data", "cpr_focal_pca_timeseries_period_1961-2017.csv")) %>%rename(PC1_small_zoo =`First Mode`,PC2_large_zoo =`Second Mode`) %>%select(-c(pca_period, taxa_used)) %>%pivot_longer(cols =starts_with("PC"), names_to ="Var", values_to ="value")```Taking PCA timeseries as proxies for those communities and evaluating them for breakpoints gives the following results.```{r}# Run breakpoints in CPR PCA# Run the regime shift testcpr_indices_rstars <- cpr_community %>%filter(year >1976) %>%mutate(EPU ="GOM",var_epu =str_c(Var, "X", EPU)) %>%split(.$var_epu) %>%map_dfr(function(.x){# detrend trend_mod <-lm(value ~ year, data = .x) .x$trend_resid <-resid(trend_mod)# cutoff length cutoff_length <-7# Get the results from that x_rstars <-rstars(data.timeseries =as.data.frame( .x[,c("year", "trend_resid")]),l.cutoff = cutoff_length,pValue =0.05,Huber =1,Endfunction = T,preWhitening = T,OLS = F,MPK = T,IP4 = F,SubsampleSize = (cutoff_length +1)/3,returnResults = T) %>%mutate(Value = .x$value,shift_direction =case_when( RSI >0~"Shift Up", RSI <0~"Shift Down",TRUE~NA)) },.id ="var_epu" ) %>%separate(var_epu, into =c("Var", "EPU"), sep ="X") %>%mutate(time =as.Date(str_c(year, "-01-01")))``````{r}# Summarise the breakpoint locationscpr_shift_points <- cpr_indices_rstars %>%filter(RSI !=0) %>% dplyr::select(time, Var, EPU, shift_direction)# Plot the breaks over the monthly dataggplot() +geom_vline(data = cpr_shift_points,aes(xintercept = time,color = shift_direction),linewidth =1.5) +geom_line(data = cpr_indices_rstars,aes(time, trend_resid),linewidth =0.4, alpha =0.5) +scale_color_gmri() +facet_grid(Var ~ EPU, labeller =label_wrap_gen(width =8), scales ="free") +scale_x_date(breaks =seq.Date(from =as.Date("1950-01-01"),to =as.Date("2020-01-01") ,by ="10 year"),limits =as.Date(c("1970-01-01", "2020-01-01")),labels = scales::label_date_short()) +theme(legend.position ="bottom",strip.text.y =element_text(angle =0)) +labs(color ="",y ="Measurement",x ="Date",title ="CPR Zooplankton Community Metrics",subtitle ="STARS Regime Shifts")```### ECOMON Community PCA Index```{r}# Abundance per 100m3 for different taxa# ecodata::zoo_regime %>% distinct(Var) %>% pull() %>% sort()ecomon_zoo <- ecodata::zoo_regimeecomon_zoo %>%filter(!str_detect(Var, "fish|clauso|gas")) %>%ggplot(aes(Time, Value)) +geom_line() +facet_grid(Var ~ EPU)# We might be able to just pull out the seven focal species and then repeat the PCA process...# Issues:# ecomon doesn't split the calanus into adult/juvenile# Para and Pseaudocalana are split# Euph also has a Euph1```### Large and Small Copepod Indexhttps://noaa-edab.github.io/tech-doc/zoo_abundance_anom.html?q=zoopl#copepod> Abundance anomalies are computed from the expected abundance on the day of sample collection. Abundance anomaly time series are constructed for Centropages typicus, Pseudocalanus spp., Calanus finmarchicus, and total zooplankton biovolume. The small-large copepod size index is computed by averaging the individual abundance anomalies of Pseudocalanus spp., Centropages hamatus, Centropages typicus, and Temora longicornis, and subtracting the abundance anomaly of Calanus finmarchicus. This index tracks the overall dominance of the small bodied copepods relative to the largest copepod in the Northeast U.S. region, Calanus finmarchicus.```{r}# This has "LgCopepods" & "SmCopepods", which could produce large/small index# ecodata::zoo_abundance_anom %>% distinct(Var) %>% pull() %>% sort()zoo_lg_small <- ecodata::zoo_abundance_anom %>%filter(Var %in%c("LgCopepods", "SmCopepods")) %>%mutate(Value =as.numeric(Value)) %>%pivot_wider(values_from ="Value", names_from ="Var") %>%mutate(small_large_index = SmCopepods - LgCopepods)ggplot(zoo_lg_small, aes(Time, small_large_index)) +geom_line() +geom_hline(yintercept =0) +facet_grid(EPU~., scales ="free") +labs(y ="Small-Large Copepod Index\n(More Large Copepods <-----> More Small Copepods)")``````{r}#| label: zooplankton, not used# # This is too many vars# ecodata::zooplankton_index %>% distinct(Var) %>% pull()# # # Not informative mechanistically# ecodata::zoo_diversity# # # BOLD move on absolute abundances# ecodata::zoo_strat_abun```### NEEDS: MCC & Lobster Predator IndicesThere are two EPU-Scale indices that we need to develop. This is the MCC index, and a lobster predator abundance index.The Gulf of Maine Coastal Current plays an important role in transporting lobster larva and their recruitment form year-to-year. The degree of "connected-ness" of the Western and Eastern portions of this current have been used in the past to inform expectations of lobster recruitment.## Local/Nearshore ShiftsConditions closer to the coast show the following long-term trends:```{r}# Temperaturefvcom_trends %>%filter(str_detect(var, "temperature"),`Trend?`,!area_id %in%c("gom_gbk", "sne")) %>%group_by(area_id) %>% gt::gt() %>% gt::tab_header(title ="FVCOM Offshore Long Term Salinity Trends",subtitle ="Monotonic Trends Evaluated by Mann-Kendall Test") %>% gt::fmt_missing(columns =everything(), missing_text ="")# Salinityfvcom_trends %>%filter(!str_detect(var, "temperature"),`Trend?`,!area_id %in%c("gom_gbk", "sne")) %>%group_by(area_id) %>% gt::gt() %>% gt::tab_header(title ="FVCOM Offshore Long Term Salinity Trends",subtitle ="Monotonic Trends Evaluated by Mann-Kendall Test") %>% gt::fmt_missing(columns =everything(), missing_text ="")```### Temperature and Salinity Breaks```{r}#| fig.height: 8# Plot the breaks over the monthly dataggplot() +# Add a vertical line using geom_segment with an arrowgeom_segment(data =filter( inshore_tempsal_shifts, str_detect(shift_direction, "Up")),linewidth =0.8,aes(x = time, xend = time, y =-Inf, yend =Inf, color = shift_direction), arrow =arrow(length =unit(0.3, "cm"), ends ="last")) +geom_segment(data =filter( inshore_tempsal_shifts, !str_detect(shift_direction, "Up")),linewidth =0.8,aes(x = time, xend = time, y =-Inf, yend =Inf, color = shift_direction), arrow =arrow(length =unit(0.3, "cm"), ends ="first")) +geom_line(data =filter(inshore_tempsal),aes(time, detrended_vals),linewidth =0.4, alpha =0.5) +scale_color_gmri() +facet_grid(area_id~Var, labeller =label_wrap_gen(width =8), scales ="free") +scale_x_date(breaks =seq.Date(from =as.Date("1950-01-01"),to =as.Date("2020-01-01") ,by ="10 year"),limits =as.Date(c("1978-01-01", "2020-01-01")),labels = scales::label_date_short()) +theme(legend.position ="bottom",strip.text.y =element_text(angle =0)) +labs(color ="",y ="Measurement",x ="Date",title ="FVCOM Inshore Scale Sal Metrics",subtitle ="STARS Regime Shifts")```#### Salinity```{r}#| fig.height: 8# Plot the breaks over the monthly dataggplot() +# Add a vertical line using geom_segment with an arrowgeom_segment(data =filter( inshore_tempsal_shifts,str_detect(Var, "Salinity"),str_detect(shift_direction, "Up")),linewidth =0.8,aes(x = time, xend = time, y =-Inf, yend =Inf, color = shift_direction), arrow =arrow(length =unit(0.3, "cm"), ends ="last")) +geom_segment(data =filter( inshore_tempsal_shifts, str_detect(Var, "Salinity"),!str_detect(shift_direction, "Up")),linewidth =0.8,aes(x = time, xend = time, y =-Inf, yend =Inf, color = shift_direction), arrow =arrow(length =unit(0.3, "cm"), ends ="first")) +geom_line(data =filter(inshore_tempsal, str_detect(Var, "Salinity")),aes(time, detrended_vals),linewidth =0.4, alpha =0.5) +scale_color_gmri() +facet_grid(area_id~Var, labeller =label_wrap_gen(width =8), scales ="free") +scale_x_date(breaks =seq.Date(from =as.Date("1950-01-01"),to =as.Date("2020-01-01") ,by ="10 year"),limits =as.Date(c("1978-01-01", "2020-01-01")),labels = scales::label_date_short()) +theme(legend.position ="bottom",strip.text.y =element_text(angle =0)) +labs(color ="",y ="Salinity Anomaly",x ="Date",title ="FVCOM Hindcast Inshore Salinity",subtitle ="STARS Regime Shifts")``````{r}inshore_tempsal_shifts %>%filter(str_detect(Var, "Salinity")) %>%group_by(Var) %>%arrange(time, area_id) %>% gt::gt() %>% gt::tab_header(title ="Inshore Salinity Regime Breaks")``````{r}# Highlight whichever regions have seen salinity regime changeggplot() +geom_sf(data =filter( study_regions, area_id %in% (dplyr::filter(tempsal_shifts, str_detect(Var, "Salinity")) %>%pull(area_id)) ),fill =gmri_cols("gmri blue"), alpha =0.4) +geom_sf(data = new_england) +geom_sf(data = canada) +coord_sf(xlim =c(-78, -66), ylim =c(35.5, 45)) +labs(title ="Salinity Regime Change Affected Areas")```#### Temperatures```{r}#| fig.height: 8# Plot the breaks over the monthly dataggplot() +# Add a vertical line using geom_segment with an arrowgeom_segment(data =filter( inshore_tempsal_shifts,!str_detect(Var, "Salinity"),str_detect(shift_direction, "Up")),linewidth =0.8,aes(x = time, xend = time, y =-Inf, yend =Inf, color = shift_direction), arrow =arrow(length =unit(0.3, "cm"), ends ="last")) +geom_segment(data =filter( inshore_tempsal_shifts, !str_detect(Var, "Salinity"),!str_detect(shift_direction, "Up")),linewidth =0.8,aes(x = time, xend = time, y =-Inf, yend =Inf, color = shift_direction), arrow =arrow(length =unit(0.3, "cm"), ends ="first")) +geom_line(data =filter(inshore_tempsal, !str_detect(Var, "Salinity")),aes(time, detrended_vals),linewidth =0.4, alpha =0.5) +scale_color_gmri() +facet_grid(area_id~Var, labeller =label_wrap_gen(width =8), scales ="free") +scale_x_date(breaks =seq.Date(from =as.Date("1950-01-01"),to =as.Date("2020-01-01") ,by ="10 year"),limits =as.Date(c("1978-01-01", "2020-01-01")),labels = scales::label_date_short()) +theme(legend.position ="bottom",strip.text.y =element_text(angle =0)) +labs(color ="",y ="Temperature Anomaly",x ="Date",title ="FVCOM Hindcast Inshore Temperature",subtitle ="STARS Regime Shifts")``````{r}ggplot() +geom_sf(data =filter( study_regions, area_id %in% (dplyr::filter(inshore_tempsal_shifts, !str_detect(Var, "Salinity")) %>%pull(area_id)) ),fill =gmri_cols("lv orange"), alpha =0.4) +geom_sf(data = new_england) +geom_sf(data = canada) +coord_sf(xlim =c(-78, -66), ylim =c(35.5, 45)) +labs(title ="Inshore Temperature Regime Shift Affected Areas")``````{r}inshore_tempsal_shifts %>%filter(!str_detect(Var, "Salinity")) %>%group_by(Var) %>%arrange(time, area_id) %>% gt::gt() %>% gt::tab_header(title ="Inshore Scale Temperature Regime Breaks")```### Days in Key Temperature RangesIn addition to breaks in absolute temperatures, there is interest in the amount of time spent in favorable (12-18C) and unfavorable conditions (20C).These use daily bottom temperatures:```{r}# Load the data for the new regionsdaily_fvcom_temps <-read_csv(str_c(lob_ecol_path, "FVCOM_processed/area_timeseries/new_regions_fvcom_temperatures_daily.csv")) %>%mutate(depth_type =if_else(area_id %in%c("gom_gbk", "SNE"), "offshore", "nearshore"),area_id =factor(area_id, levels = areas_northsouth) )# # Monthly Salinity# monthly_fvcom_sal <- read_csv(# str_c(lob_ecol_path, "FVCOM_processed/area_timeseries/new_regions_fvcom_salinity_monthly_gom3.csv")) %>% # mutate(# depth_type = if_else(area_id %in% c("gom_gbk", "SNE"), "offshore", "nearshore"),# area_id = factor(area_id, levels = areas_northsouth)# )# The degday functions were built with the ability to model daily cycles# the sine and triangle methods accomodate this# since we only have daily data the simple average is probably the way to do itthresh_low <-10thresh_up <-18# library(degday)# Get Monthly Totals in Rangesdd_monthly <- daily_fvcom_temps %>%mutate(year = lubridate::year(time),month = lubridate::month(time),opt_btemp =if_else(between(bottom_t, 10,18), 1, 0),stress_btemp =if_else(bottom_t >18, 1, NA),cold_btemp =if_else(bottom_t <10, 1, NA)) %>%group_by(area_id, year, month) %>%summarise(across(ends_with("temp"), ~sum(.x, na.rm = T)),.groups ="drop") %>%pivot_longer(cols =ends_with("temp"), names_to ="var", values_to ="totals") %>%mutate(time =as.Date(str_c(year,"-01-01")) +months(month-1),area_id =factor(area_id, levels = areas_northsouth),var =case_match( var,"opt_btemp"~"Preferred Bottom Temperatures 10-18C", "stress_btemp"~"Heat Stress Conditions >18C","cold_btemp"~"Below Preferred Conditions <10C") )``````{r}# Do annual, Monthly values looked insanedd_annual <- dd_monthly %>%group_by(year, area_id, var) %>%summarise(across(totals, sum)) %>%mutate(var_area =str_c(var, area_id, sep ="X"),time =as.Date(str_c(year, "-01-01")))# Plot themdd_annual %>%ggplot() +geom_area(aes(year, y = totals, fill = var)) +scale_fill_manual(values =c("lightblue", "#ea4f12", "#057872")) +facet_grid(area_id~var) +scale_x_continuous(expand =expansion(add =c(0,0))) +theme(strip.text.y =element_text(angle =0),legend.position ="bottom") +guides(fill =guide_legend(nrow =2,title.position ="top",title.hjust =0.5))+labs(y ="Days in Range",fill ="Daily Temperature Conditions", color ="",title ="FVCOM Bottom Temperature Degree-Days")``````{r}# # Remove monthly averages, and trendsdd_annual_detrended <- dd_annual %>%split(.$var_area) %>%map_dfr(function(.x){# Detrend .x <- .x %>%arrange(year) %>%mutate(yr_num =as.numeric(year))# annual trend + monthly average trend_mod <-lm(totals ~ yr_num, data = .x)# save the results .x <- broom::augment(x = trend_mod) %>%rename(trend_fit = .fitted,trend_resid = .resid) %>%full_join(.x, join_by(yr_num, totals)) %>%mutate(trend_resid =if_else(is.na(totals), NA, trend_resid))return(.x)}) %>%mutate(time =as.Date(str_c(year, "-01-01")))# Plotdd_annual_detrended %>%filter(area_id %in%tolower(c("GOM_GBK", "SNE"))) %>%ggplot(aes(time, trend_resid, fill = var, color = var)) +geom_col() +geom_smooth(method ="loess", linewidth =0.6) +scale_color_manual(values =c("lightblue", "#ea4f12", "#057872")) +scale_fill_manual(values =c("lightblue", "#ea4f12", "#057872")) +facet_grid(area_id~var, scales ="free", labeller =label_wrap_gen(width=8)) +guides(color =guide_legend(nrow =3)) +labs(title ="EPU Scale Ocean/Climate Metrics - Detrended",y ="Departure from Long-Term Trend (days)")``````{r}# Run the regime shift testtemp_range_rstars <- dd_annual_detrended %>%# temp_range_rstars <- dd_annual %>% split(.$var_area) %>%map_dfr(function(.x){# Seven years cutoff_length <-7# Get the results from that x_rstars <-rstars(data.timeseries =as.data.frame( .x[,c("time", "trend_resid")]),l.cutoff = cutoff_length,pValue =0.05,Huber =1,Endfunction = T,preWhitening = T,OLS = F,MPK = T,IP4 = F,SubsampleSize = (cutoff_length +1)/3,returnResults = T) %>%mutate(Value = .x$totals,shift_direction =case_when( RSI >0~"Shift Up", RSI <0~"Shift Down",TRUE~NA)) },.id ="Var" )%>%separate(Var, into =c("Var", "area_id"), sep ="X") %>%mutate(area_id =factor(area_id, levels = areas_northsouth))```### Temperature Suitability ShiftsThe results can be seen below:```{r}#| fig.height: 8# Summarise the breakpoint locationstemp_suit_shift_points <- temp_range_rstars %>%filter(RSI !=0) %>% dplyr::select(time, Var, area_id, shift_direction)# Plot the breaks over the monthly dataggplot() +# Add a vertical line using geom_segment with an arrowgeom_segment(data =filter( temp_suit_shift_points, str_detect(shift_direction, "Up")),linewidth =0.8,aes(x = time, xend = time, y =-Inf, , yend =Inf, color = shift_direction), arrow =arrow(length =unit(0.3, "cm"), ends ="last")) +geom_segment(data =filter( temp_suit_shift_points, !str_detect(shift_direction, "Up")),linewidth =0.8,aes(x = time, xend = time, y =-Inf, , yend =Inf, color = shift_direction), arrow =arrow(length =unit(0.3, "cm"), ends ="first")) +geom_line(data = temp_range_rstars,aes(time, Value),linewidth =0.4, alpha =0.5) +scale_color_gmri() +facet_grid(area_id ~ Var, labeller =label_wrap_gen(width =8), scales ="free") +scale_x_date(breaks =seq.Date(from =as.Date("1950-01-01"),to =as.Date("2020-01-01") ,by ="10 year"),limits =as.Date(c("1970-01-01", "2020-01-01")),labels = scales::label_date_short()) +theme(legend.position ="bottom",strip.text.y =element_text(angle =0)) +labs(color ="",y ="Measurement",x ="Date",title ="Lobster Thermal Preferences - Detrended",subtitle ="STARS Regime Shifts")```Based on the annual totals, we see limited breakpoints in suitable thermal habitat.```{r}temp_suit_shift_points %>%group_by(Var) %>%arrange(time, area_id) %>% gt::gt() %>% gt::tab_header(title ="Temperature Suitability Scale Breaks")```And restricted to these areas```{r}ggplot() +geom_sf(data =filter( study_regions, area_id %in% temp_suit_shift_points$area_id),fill =gmri_cols("lv orange"), alpha =0.4) +geom_sf(data = new_england) +geom_sf(data = canada) +coord_sf(xlim =c(-78, -66), ylim =c(35.5, 45)) +labs(title ="Affected Areas")```## Summary Figures / TablesDo the trend evaluation for everything:Summary table should have the region, the variable, whether there is a trend or not, what the rate is, then whether there have been breakpoints, and when they were (mm/yyyy).